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Abstract 

We  develop  a  diffuse  interface  model  to  investigate  the  three  dimensional  coarsening  in  thin 
films.  In  this  model,  both  strong  surface  anisotropy  with  linear  Willmore  regularization, 
elastic  stresses  and  deposition  are  included.  The  governing  equation  for  the  phase  field 
parameter  is  a  sixth  order  Cahn-Hilliard  Equation  due  to  the  presence  of  surface  anisotropy 
and  the  linear  Willmore  regularization.  The  simulated  system  is  assumed  to  be  in  mechanical 
equilibrium  with  misfit  in  the  film  generated  by  lattice  mismatch  in  the  substrate.  Thus, 
the  Cauchy-Navier  equations  are  solved  for  elastic  displacements  which  lead  to  the  elastic 
energy.  Both  the  Cahn-Hilliard  equation  and  the  Cauchy-Navier  equations  are  solved  with 
an  non-stiff,  adaptive  nonlinear  multigrid  method.  Simulation  results  of  coarsening  in  three 
dimensions  with  different  strengths  of  the  surface  anisotropy,  misfit  strain,  and  deposition 
rates  are  shown.  Comparison  and  analysis  of  these  results  help  to  explain  their  influence  on 
coarsening  processes  in  thin  films. 


Introduction 

It  is  known  that  both  surface  anisotropy  and  elastic  strains  have  an  important  influence  on 
the  morphological  evolution  of  thin  films.  The  effects  of  strong  surface  anisotropy  have  been 
studied  in  great  detail  in  the  problem  of  crystal  growth [1].  In  the  case  of  an  isotropic  surface 
energy,  a  crystal  usually  shows  a  spherical  shape  with  all  crystal  orientations  present.  In  con¬ 
trast,  strong  surface  anisotropies  result  in  missing  crystal  orientations  and  enable  formation 
of  facets  at  the  surface [2].  However,  non-convexity  of  this  strong  surface  anisotropy  energy 
leads  to  ill-posedness  of  the  Cahn-Hilliard  equation [3]  [4].  Hence,  Willmore  regularization  has 
been  proposed  to  resolve  this  issue  by  smoothing  out  the  transitions  near  the  sharp  corners 
at  the  corrugated  surface  [5]. 

Besides  surface  anisotropy,  elastic  strains  also  play  an  important  role  in  the  morphological 
evolution  of  thin  films.  Due  to  the  difference  in  crystal  structures  between  the  substrate  and 
film,  misfit  strains  arise.  The  elastic  energy  stored  in  the  thin  films  owing  to  these  misfit 
strains  can  be  relieved  by  developing  pits,  islands  and  quantum  fortresses [6].  Formation  of 
these  structures  increase  the  free  surface  area  to  aid  the  relaxation  of  elastic  energy,  though 
at  the  cost  of  increasing  the  surface  free  energy [7].  Competition  between  the  elastic  energy 
and  the  surface  free  energy,  combined  with  deposition  of  materials  at  the  surface  dominate 
the  process  of  morphological  evolution  of  thin  films. 

In  this  paper,  a  diffuse  interface  model  is  used  to  study  the  coarsening  of  thin  films  in  three 
dimensions  with  different  strengths  of  the  surface  anisotropy,  misfit  strains,  and  deposition 
rates.  Formulation  of  governing  equations  is  presented  in  section  2.  Numerical  methods 
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are  given  in  section  3.  Preliminary  results  with  different  strengths  in  surface  anisotropy  are 
presented  and  discussed  in  section  4. 


Formulation 

The  thermodynamic  system  considered  here  consists  of  a  thin  film  phase  and  a  vapor  phase. 
A  phase  variable  c  is  introduced  such  that  c  =  0.0  is  the  vapor  phase  and  c  =  1.0  is  the  thin 
film  phase.  Both  strong  surface  anisotropy  and  linear  Willmore  regularization  have  been 
introduced  to  the  model  to  give  the  surface  a  corrugated  morphology  with  smooth  corners. 
Elastic  stresses  arise  from  the  misfit  strains  are  also  included. 

Thus  the  total  free  energy  of  this  system  is  given  by: 

E  =  /  i7(n)(/(c)  +  ^|Vc|2)dX  (1) 

Jn  e  2 

+ 

+  [  W(c ,  E)dX, 

Jn 

where 

/(c)  =  -^j4c2(1  -  C)2.  (2) 

/(c)  is  a  double  well  chemical  free  energy  density  with  c  =  0.0  and  c  =  1.0  being  the  two 
minima,  and  A  is  a  parameter  related  to  the  height  of  the  energy  barrier;  ~  |Vc|2  is  the 
gradient  energy  density  which  takes  into  account  the  non-local  interactions  due  to  the  phase 
field  variation  across  the  interface; 


7(n)  =  l  +  al/n), 

(3) 

r(n)  =  4X>4-3, 

d  =  2,3 

(4) 

i—1 

Vc 

n  ~  1W 

(5) 

where  7(n)  is  the  function  of  surface  anisotropy  with  n  being  the  normal  vector  at  the 
interface;  a  is  dimensionless  and  defines  the  strength  of  the  anisotropy,  and  r(n)  as  given 
in  Equation  (2)  indicates  that  the  anisotropy  is  four-fold  symmetric  with  n*  being  each 
component  of  the  normal  vector  n.  As  shown  in  equation  (1),  the  anisotropy  surface  energy 
7(11)  times  both  the  chemical  free  energy  denstiy  and  the  gradient  energy  density.  It  has 
been  shown  in  [8]  that  the  thickness  of  interfaces  from  this  new  formulation  is  uniform  and 
independent  of  the  crystal  orientation. 

The  second  term  given  in  Equation  (1)  is  the  linear  Willmore  regularization  term  and  (3 
is  regularization  parameter.  The  sharp  corners  due  to  formation  of  facets  will  be  smoothed 
out  due  to  the  higher-order  derivative  in  this  term.  The  last  term  W (c,  E)  is  the  elastic 
energy  given  by, 

I'V. /••"!  1  £  Wj-Mc))-  (6) 

^  i,j= 1 
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Accordingly,  we  define 


d 

Tij  =  E  Cm{c){Ekl-Skle{c)),  (7) 

i,j= 1 

Rkl  =  “I-  (8) 

e(c)  =  rjq(cy,q(c)  =  c2(3-  2c),  (9) 

C««  =  C0m+q{c){CljM-C°m),  (10) 


where  and  are  linear  stress  and  strain  tensor  related  by  Hooke’s  Law  as  given  in 
Equation  (7);  u  are  elastic  displacements;  e(c)  is  the  misfit  strain  arising  from  the  mismatch  of 
crystal  lattice  with  the  magnitude  given  by  77;  the  cubic  elastic  stiffness  tensor  Cijki  is  defined 
to  be  an  interpolation  of  the  stiffness  tensors,  C^kl  and  Cfjkl,  of  the  two  phases  weighted  by 
q(c).  For  the  vapor  phase,  C^kl  are  assigned  with  values  three  orders  of  magnitude  lower 
than  those  of  the  film  phase. 

Evolution  of  the  phase  field  is  determined  by  the  sixth-order  regularized,  anisotropic 
Cahn-Hilliard  equation  as  shown  below, 

V  •  (M(c)Vju)  +  S(c,  n),  (11) 

7(n)/'(c)  -  e2V  •  (Vc  +  VA(n)Vc)  -  0Av  +  W'(c,  E),  (12) 

e2  Ac,  (13) 

where  the  generalized  chemical  potential  /j,  is  the  variational  derivative  of  the  total  free 
energy  density;  M(c)  is  the  surface  mobility  given  by  Mqc(1  —  c)  with  Mq  being  the  mobility 
coefficient;  S(c,  n)  =  VdQRc2(l  —  cyn^  is  the  deposition  term  localized  at  the  film- vapor 
interface[9]  with  Vd  being  the  deposition  rate,  Q  is  a  scale  factor  to  match  the  sharp  interface 
result,  and  R  is  a  random  number(0.9  <  R  <  1.1).  For  computational  convenience,  the 
Cahn-Hilliard  equation  is  split  into  three  second-order  equations  with  c,  ji  and  v  being  the 
three  dependent  variables. 

Elastically,  the  film  phase  and  the  vapor  phase  are  assumed  to  be  in  mechanical  equilib¬ 
rium.  The  conditions  of  mechanical  equilibrium  are  used  to  determine  the  elastic  fields. 

d 

y  o,  d  =  2,3.  (14) 

3= 1 

This  quasi- stationary  assumption  combined  with  stresses  and  strains  given  in  Equation  (7) 
lead  to  the  Carchy-Navier  equations  for  elastic  displacements.  Elastic  energy  is  determined 
via  Equation  (6)  once  the  displacements  are  found  from  these  Carchy-Navier  equations. 

For  these  governing  equations,  periodic  boundary  conditions  are  used  in  the  horizontal 
directions  for  each  dependent  variable.  Neumann  boundary  conditions  are  used  for  the  Cahn- 
Hillard  equations  in  the  vertical  direction  of  the  computational  cell;  Traction-free  boundary 
conditions  and  Dirichlet  boundary  conditions  are  used  for  the  displacements  at  the  top  and 
bottom  of  the  computational  cell,  respectively. 

N  ondimensionalization 

All  governing  equations  are  nondimensionalized  for  computational  convenience.  The  scaled 
position  and  scaled  time  are  defined  to  be  x  =  ^  and  t  =  ^  where  the  characteristic  Length 


8c 

dr 

[i  = 
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L  is  chosen  to  be  20 ran  and  the  characteristic  time  T  is  given  as  below: 


T  = 


L2 

M0A’ 


(15) 


The  chemical  free  energy  density  is  normalized  by  the  parameter  A,  and  M(c)  by  Mq. 
The  stiffness  coefficients  are  normalized  by  C44.  The  misfit  strain  and  the  elastic  energy 
density  are  normalized  by  77  and  772C|4,  which  are  usually  used  as  the  characteristic  scales  for 
the  elastic  energy [9].  Other  non-dimensional  quantities  characterizing  this  physical  system 
are  listed  below, 


e  .5  VdL2  7_£CL 

LVA I?'Vd  M0A ’  A 


(16) 


where  Z  is  a  parameter  relating  to  the  strength  of  the  elastic  effects.  Most  of  the  nondi- 
mensionalized  quantities  are  denoted  with  an  over  bar  correspondingly.  The  non-dimensional 
governing  equations  can  be  derived  by  equating  the  operators  and  sources  shown  by  equations 
(20)- (24)  and  (25)- (29)  in  the  next  subsection,  respectively. 

The  physical  values  for  the  characteristic  length  and  time,  and  the  calculated  values  for 
the  normalized  quantities  except  those  for  the  deposition  term  are  listed  in  Table  1. 


paramter 

L 

T  e 

P 

z 

value 

20(nm) 

10-100(s)  2.4  x  ltr2 

0.01 

0.125-0.25 

Table  1:  Physical  values  and  calculated  values  used  in  simulation 

Numerical  methods 

The  governing  equations  given  in  equations  (11)  and  (14)  are  solved  by  adaptive  nonlinear 
multigrid  methods  developed  by  Wise,  Kim  and  Lowengrub[4].  Briefly,  the  sixth  order  Cahn- 
Hilliard  equation  with  linear  Willmore  regularization  is  split  into  three  second  order  equations 
as  shown  in  equation  (11).  These  three  equations  and  the  elastic  Cauchy-Navier  equations 
are  solved  simultaneously  using  an  adaptive  multilevel  multigrid  method  to  determine  the 
profile  of  the  composition  field  c  and  the  displacement  field  u.  Second  order  center  finite 
differences  are  used  for  the  spatial  discretization  and  a  gradient  stable  scheme  is  adopted  for 
the  time  discretization. 

For  the  gradient  stable  scheme,  the  normalized  chemical  free  energy  density  is  split  into 
a  contractive  term  and  an  expansive  term  as  shown  below, 


= 

tic)  -  tic), 

(17) 

tic) 

=  -(c-  , 

2  64 

(18) 

tic) 

8'  2 

(19) 

A  full  approximation  scheme  (FAS)  multigrid  method  is  applied  to  solve  the  nonlinear 
governing  equations  and  self  adaptive  block-structured  Cartesian  mesh  refinement  is  used  to 
resolve  narrow  interfacial  layers.  For  the  FAS  method,  equations  are  split  into  operator  (N) 
and  source  (S)  terms [4]  with  the  components  of  the  operator  being 

Afi  =  cn+1  -  dt\7  •  M(cn+1’fc)V/T+1,  (20) 

N2  =  pn+1  -  fc{cn+1)  +  nun+1  -  (3Ailn+\  (21) 
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and  those  of  source  terms  being 


51  =  cn  +  VdQR{ef(l-cnfnz, 

52  =  —f'e(cn)  +  otT(n)f'(cn)  —  e2V  ■  J4(n)Vc™+1  +  W' (cn+1 ,  En+1) , 

53  =  0, 

54  =  ^-[(C'11  +  C'12)e(c"+1)], 

55  =  [(C22  +  C'12)e(cn+1)j  , 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


The  operator  and  source  terms  for  the  Cauchy-Navier  equations  are  only  shown  in  the  2D 
case.  The  3D  case  is  analogously  defined.  More  details  about  the  numerical  method  can  be 
found  in  the  work  of  Wise  et  al  [4] . 

The  size  of  the  computational  cell  is  set  to  be  6.4  x  6.4  x  3.2.  Two  refinement  levels  are 
used  and  the  root  level  is  a  grid  of  size  64  x  64  x  32.  Thus  the  finest  mesh  size  is  h  =  0.025. 
The  time  step  used  in  the  simulation  is  1.6  x  10-2.  The  tolerance  set  for  the  convergence  of 
the  multigrid  solver  is  1.0  x  10-6. 

Compared  with  the  Crank-Nicholson  method  for  time  discretization,  the  gradient  stable 
method  can  allow  a  time  step  which  is  16  times  of  that  in  the  former  case.  In  the  3D 
simulations  presented  below,  the  estimated  execution  time  is  successfully  reduced  from  8-10 
months  to  about  one  month. 


Preliminary  Simulation  Results 


In  this  section,  preliminary  simulation  results  are  shown  that  describe  the  influence  of  the 
strength  of  surface  anisotropy  without  material  deposition.  Fig.  1  shows  the  morphological 
evolution  of  thin  films  with  an  initial  sinusoidal  perturbation  which  is  given  by, 


c(xltx2,x3) 


1  1.0  —  tanhfe  —  r(x  i,  x2)) 

2  2 Til 


(30) 


where 

r(xi,  X2)  =  1.6  +  0.06szn(^^1)  sin(^J^2).  (31) 

6.4  6.4 

The  magnitude  of  this  sinusoidal  perturbation  is  0.06.  The  surface  anisotropy  parameter  a 
is  set  to  be  a  =  0.2  in  Fig.  1(a)  and  (b),  and  a  =  0.3  in  (c)  and  (d),  respectively.  However, 
in  both  cases  the  strengths  of  the  elasticity  are  the  same  with  Z  =  0.25.  A  comparison  of 
Fig.  1(a)  and  (c),  and  (b)  and  (d)  at  the  non-dimensional  times  t  =  1.6  and  t  =  6.4,  show 
that  the  pyramids  formed  in  (c)  and  (d  ) present  sharper  corners  and  edges  when  the  surface 
anisotropy  is  stronger.  In  Fig.  1(c),  pits  are  also  found  to  form  at  the  top  of  the  pyramids 
to  help  decrease  the  surface  anisotropy  energy. 

Fig.  2  shows  the  morphological  evolution  of  thin  films  with  an  initial  random  perturbation 
as  given  by  equation  (30)  with  r(xi,x2)  replaced  by, 


r(xi,x2)  =  1.6  +  0.0001£(£i,  x2). 


(32) 
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(a) 


(b) 


(c)  (d) 

Figure  1:  Morphological  evolution  of  thin  films  with  an  initial  sinusoidal  perturbation,  and 
Z  =  0.25:  (a)  a  =  0.2,  t  =  1.6;  (b)  a  =  0.2,  t  =  6.4;  (c)  a  =  0.3,  t  =  1.6;  (d)  a  =  0.3, 
t  =  6.4. 

where  C{xux2)  is  a  random  number  and  the  magnitude  of  the  perturbation  is  10-4.  Other 
physical  parameters  are  set  to  be  the  same  as  in  the  former  case.  Compared  with  those 
shown  in  Fig.  2(a)  and  (b),  the  morphological  evolution  in  Fig.  2(c)  and  (d)  proceeds  much 
faster  due  to  the  stronger  surface  anisotropy.  At  time  t  =  1.6,  the  surface  in  Fig.  2(a)  with 
a  =  0.2  is  still  close  to  being  planar,  while  that  in  Fig.  2(c)  has  already  become  highly 
corrugated.  The  difference  shown  by  these  preliminary  results  suggests  that  at  the  early 
times,  surface  anisotropy  is  the  dominant  driving  force  for  morphological  evolution,  since  a 
corrugated  surface  can  greatly  lower  the  surface  anisotropy  energy  initially. 

In  Fig.  3,  longer-time  simulation  results  with  a  =  0.2  and  Z  =  0.25  are  shown  for  both 
initial  perturbations.  Fig.  3(a)  and  (b)  demonstrate  well- developed  domes  and  valleys  for 
the  surface  of  the  thin  film  with  a  sinusoidal  perturbation.  As  time  evolves,  the  valleys 
are  supposed  to  deepen  to  aid  the  relaxation  of  elastic  stresses  and  strains,  thus  lower  the 
total  elastic  energy.  No  coarsening  seems  to  occur  in  these  two  plots  since  the  domes  are 
approximately  the  same  size  due  to  the  initial  perturbation.  However,  in  Fig.  3(c)  and  (d), 
coarsening  of  the  domes  is  observed.  This  process  lowers  the  total  elastic  energy.  Due  to 
the  use  of  the  gradient  stable  method,  the  total  execution  time  for  these  two  longer-time 
simulations  is  about  2  weeks. 
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(a) 


(b) 


(c)  (d) 

Figure  2:  Morphological  evolution  of  thin  films  with  an  initial  random  perturbation,  and 
Z  =  0.25:  (a)  a  =  0.2,  t  =  1.6;  (b)  a  =  0.2,  t  =  6.4;  (c)  a  =  0.3,  t  =  1.6;  (d)  a  =  0.3, 
t  =  6.4. 

Summary 

In  this  paper,  a  diffuse  interface  model,  including  strong  surface  anisotropy  with  linear  Will- 
more  regularization,  elastic  stresses  and  deposition,  is  developed  to  simulate  the  three  dimen¬ 
sional  coarsening  in  thin  films.  The  high-order  governing  equations,  i.e,  the  Cahn-Hilliard 
equation  and  the  Cauchy-Navier  equations  are  solved  with  an  non- stiff,  adaptive  nonlinear 
multigrid  method.  A  gradient  stable  method  is  applied  for  the  time  discretization  which  sig¬ 
nificantly  reduces  the  total  execution  time  of  the  three  dimensional  simulations.  Preliminary 
simulation  results  of  coarsening  in  three  dimensions  with  different  strengths  of  the  surface 
anisotropy  suggests  that  surface  anisotropy  is  the  dominant  driving  force  for  morphological 
evolution  initially  since  a  corrugated  surface  can  greatly  lower  the  surface  anisotropy  energy. 
More  simulation  results  will  be  provided  in  the  future  to  further  characterize  the  coarsening 
processes  in  thin  films. 
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(a) 


(b) 


(c)  (d) 

Figure  3:  Longer-time  morphological  evolution  of  thin  films  with  a  =  0.2,  Z  =  0.25:  and 
with  an  initial  sinusoidal  perturbation  (a)  t  =  11.2,  (b)  t  =  16;  and  with  an  initial  random 
perturbation  (c)  t  —  11.2,  (d)  t  =  16. 
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